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ABSTRACT 

We examine the formation of groups of multiple supermassive black holes (SMBHs) in gas- 
poor galactic nuclei due to the high merger rate of galaxies at high redshifts. We calculate 
the relative likelihood of binary, triple, and quadruple SMBH systems, by considering the 
timescales for relevant processes and combining merger trees with N-body simulations for the 
dynamics of stars and SMBHs in galactic nuclei. Typical haloes today with mass Mo ~ 10'^ 
Mq have an average mass M~=h = 5 x lO" M© at z ~ 6, while rare haloes with current mass 
Md > 10'^ Mq have an average mass M,=6 = 5 x 10^^ M© at that redshift. These cluster-size 
haloes are expected to host single galaxies at z ~ 6. We expect about 30% galaxies within 
haloes with present-day mass Mo ~ 10'"* Mq to contain more than two SMBHs at redshifts 
2 < z < 6. For larger present-day haloes, with Mo > 10'^ Mq, this fraction is almost 60%. The 
existence of multiple SMBHs at high redshifts can potentially explain the mass deficiencies 
observed in the cores of massive elliptical galaxies, which are up to 5 times the mass of their 
central BHs. Multiple SMBHs would also lead to an enhanced rate of tidal disruption of stars, 
modified gravitational wave signals compared to isolated BH binaries, and slingshot ejection 
of SMBHs from galaxies at high speeds in excess of 2000 km s ' . 

Key words: black hole physics - galaxies: nuclei - galaxies:evolution - galaxies:high redshift 
- galaxies:kinematics and dynamics 



1 INTRODUCTION 

Most local gala xies host super massive bla ck holes (SMBHs ) 
at their centres jRichstone et al.l 1 1998 . Ferrarese & Ford 120050 . 
The SMBH mass Mbh is correlated with properties of the 
spheroidal nucleus of the host galaxy, such as velocity disper- 
sion ( Ferrarese & Merritt 2000l: iGebhardt et alj 120 00: Ferrarese 
'2002; Gultekin et al. 20091 ) and lumin osity (Magorrian et al. 1998;; 
McLure & Dunlog .2002 ; iMarconi & Hunt .2003. ; .Gultekin et all 
2009[). Detection of bright quasars at redshifts z > 6 JFan et al.l 
200 ll ; iMortlock et aLlbOl lt) suggests that SMBHs with masses as 
high as ~ 2 X 10' Mq already existed at z ~ 7. In the stan- 
dard ACDM cosmological model, growth of galaxies is hierar- 
chical and galaxy mergers are expected to be particularly fre- 
quent at redshifts z ~ 6-20. As galaxies merge, their central 
SMBHs can grow through coalescence and accretion of gas. It 
is commonly postulated that SMBHs at lower reds hifts grew out 
of se e d black holes (BH s ) in t h e first galaxies JLoeb & Rasid 
I994I; lEisenstein & Loebl Il995l; iKauffmann & Haehneltl | 200(]| ; 



Menou et alj I2OOII; iBromm & Loed |2003|; IVolonteri et al.1 I2OO3I; 



that any binary black hole system, which inevitably forms in a 
galaxy's merger history, coalesces on a short time-scale. How- 
ever, the evolution of SMBH binaries is a complex open prob- 
lem and it is unclear if a binary can merge within a Hubble time 
JMerritt & Milosavlievic 2005). One expects that during a merger 
event of two galaxies, the dyna mics of their constitu ents SMBHs 
would proceed in three stages ffiegelman et al.ll 19801) . In the first 
stage, the SMBHs sink to the centre of the gravitational potential 
of the merger remnant by dynamical friction and form a gravita- 
tionally bound binary. The newly-formed binary continues to lose 
energy and angular momentum through its global gravitational in- 
teraction with many stars until the separation between the SMBHs 
reduces to a value at which the dominant mechanism of energy loss 
is the 3-body interaction between the binary and individual stars. 
This is the second stage of the binary's evolution, and is known 
as the 'hard stage.' The precise definition of a hard SMBH binary 
varies in the literature, but it is commonly assumed that the binary 
becomes hard when its semi-major axis a reaches a value given by 
( IYJ2OO2I) 



Hopkins et al.ll200d ; lTanaka & Haimadl20 09) 

Existing merger tree models are based on the assumption 
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where stars in the galactic nucleus are assumed to have a one- 
dimensional velocity dispersion a, and m denotes the mass of the 
lighter SMBH. Finally, once the SMBH separation decreases to a 
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small-enough value, gravitational wave emission becomes the dom- 
inant mode of energy loss and the SMBHs coalesce rapidly. This is 
the third stage of the SMBH binary evolution. The value of semi- 
major axis a at which the coalescence ti me scale due to gravita - 
tional wave emission alone is t is given by jPetersI 1 964l ; lLoebl20 101) 



a(0 ^fl,„(0 = 4.3 X 10-3 (^j ( 



M 



2 X lO^Mcr 



3/4 



PC 
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where M is the total mass of the binary, and we have considered two 
SMBHs with mass 10* Mq each on a circular orbit (with shorter 
time scale at increasing eccentricity). Gravitational wave emission 
takes over as the dominant mode of energy loss when a = a,„(ti,), 
where /;, is the hardening time scale. 

Among these three stages of evolution of an SMBH binary, 
the largest uncertainty in the binary's lifetime originates from the 
hard stage, which can be the slowest of the three stages since 
the binary quickly ejects all low angular momentum stars in its 
vicinity, thus cutting off it s supply of stars. This is know n as 
the "final pars ec pr oblem" jMilosavlievic & Merrittll2003a) . For 
example, \Y^ ( l2002f) studied coalescenc e of SMBH binaries in 
a sample of galaxies observed by Fabe r et al.l ( 119971) and found 
that spherical, axisymmetric or weakly triaxial galaxies can all 
have long-lived binary SMBHs that fail to coalesce. Similarly, 
iMerritt & Milosavlievig ( 120051) found that the time spent by a bi- 
nary is less than 10'" yr only for binaries with very low mass ratios 
(< 10"^)I3 Furthermore, Merritt & Milo savlievia (|2005) showed 
that a binary may not be able to interact with all the stars in its loss 
cone, thereby increasing the time spent in the hard stage even fur- 
ther; they found that in a nucleus with a singular isothermal sphere 
stellar density profile, an equal-mass binary will stall at a separation 
of a a; ai,/2.5, where we have defined a;, in Equation (|T](. The final 
separation is expected to be even higher for galaxies with shallower 
density profiles. 

Several ways have been discussed in the literature to effi- 
ciently extract energy and angular momentum from a hard SMBH 
binary a nd overcome the final parsec problem. An example is the 
work bv lArmitage & NataraianI ( 120020 ■ who suggested that gas can 
catalyse the coalescence of a hard SMBH binary by serving as 
an effective sink for the binary's angular momentum. In particu- 
lar, they found that a binary with a separation of 0. 1 pc embed- 
ded in a gaseous accretion disk would merge in 10^ years with- 
out significant enhancem ent in the gas accretion rate. Similarly, 
lEscala et all J2004 12005|) found that in SPH simulations, clouds 
of hot gas (Jgas ~ Tyijiai) can induce decay of orbits of embedded 
binary point masses due to gravitational drag. A caveat to these 
studies is that feedback from gas accretion onto the SMBHs can re- 
move the rest of the gas from the merger remnant before the binary 
coalesces. However, stellar dynamical processes could als o acceler- 
ate bin ary coalescence, without gas. For example, Merri tt & Poonl 
( |2004[) considered the effect of chaotic orbits in steep triaxial po- 
tentials. They found that stars are supplied to the central black hole 
at a rate proportional to the fifth power of the stellar velocity dis- 
persion and that the decay rate of a central black hole binary would 
be enhanced even if only a few percent of the stars are on chaotic 
orbits, thus solving the final parsec problem. As another example, 
it was suggested that a third SMBH closely interacting with a hard 



' However, for sucli low mass ratios the time taken by the lighter black hole 
to reach the galactic nucleus due to dynamical friction is itself expected to 
exceed the Hubble time. 



SMBH binary can reduce the binary separation to a small value ei- 
ther due to the eccentricit y oscillations indu ced in the binary via the 
Kozai-Lidov mechanism dBlaes et al.ll2002n or due to repopulation 
of the binary's loss cone due to the pertu rbation in the large-scal e 
potential ca used by the third black hole (Hoff man & Loebll2007n . 
iBlaes et alj (_2PP2) found that the merger time scale of an inner cir- 
cular binary can be shortened by as much as an order of magnitude, 
and that general relativistic precession does not destroy the Kozai- 
Lidov effect for hierarchical triples that are compact enough. 

In summary, there is substantial uncertainty in the current 
understanding of the evolution of binary SMBHs. Clearly, if the 
SMBH binary coalescence time is longer than the typical time be- 
tween successive major mergers of the galaxy, then more than two 
SMBHs may exist in the nucleus of a merger remnant. We study 
this possibility in this paper. We calculate the relative likelihood 
of binary, triple, and quadruple SMBH systems, by considering 
the timescales for relevant processes and combining galaxy merger 
trees with direct-summation N-body simulations for the dynamics 
of stars and SMBHs in galactic nuclei. An obvious question regard- 
ing galactic nuclei with multiple SMBHs is whether such systems 
can be long-lived. We consider this question here. Finally, systems 
with multiple SMBHs are likely to be interesting because of obser- 
vational effects involving their effect on the properties of the host 
bulge, the enhancement in the rate of tidal disruption of stars, their 
associated gravitational wave and electromagnetic signals, and the 
slingshot ejection of SMBHs at high speeds. We study some of 
these effects. 

In ^we review previous results on galactic nuclei with more 
than two SMBHs. We present simple analytical arguments regard- 
ing the formation and evolution of such systems in ^and|4] Details 
about our numerical simulations are described in ^ with their re- 
sults shown in ^ We consider the observational signatures of our 
findings in ^ Finally, we discuss and summarise our primary find- 
ings in ^ 



2 PREVIOUS WORK 

Galactic nuclei w ith multiple SMBHs were first studied by 
ISaslaw et al.l J1974) . who computed orbits of three and four SMBH 
systems by sampling the parameter space of the problem. They 
showed that if an infalling SMBH is lighter than the components of 
the pre-existing binary, then the most probable outcome is a sling- 
shot ejection in which the infalling SMBH escapes at a velocity 
that is about a third of the orbital velocity of the binary. IValtonenl 
( Il976h further showed that the ejection velocity can be significantly 
enhanced if drag forces due to gravitational radiation are accounted 
for in the three-body dynamics. The formation of systems with 
three or four SMBHs in a hierarchical merger of smooth galac- 
tic potentials was first studied by Mikkola & Valtonen ( 1990) and 
IValtonen et alj ( ll994l) with the objective of understanding the struc- 
ture of extragalactic radio sources. This li ne of work was exte nded 
to binary-binary scat t ering of SMBHs by iHeinamakil <200ll) , and 
by iHoffman & Loebl ( 120071) ■ who studied repeated triple interac- 
tions in galactic nuclei. Both of these studies used cosmologically 
consistent initial conditions based on the extended Press-Schechter 
theory . Systems w i th a l arger number of black hol es were stud- 
ied bv lHut&Reei ( Il992h and lXu & Ostriked d 19941) using simple 
analytical models and numer ical calcu lations of massive particles 
in smooth galactic potentials. IXu & Ostriker ( 1994) concluded that 
the most-likely outcome in these cases is one in which most black 
holes are ejected and the galactic center is left with zero, or one. 
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or two black holes. Finally, full N-body simulations of galactic 
nuclei with constituent SM BHs were performed for th e case of 
two su cces sive mergers by iMakino & Ebisuzakil ( Il996r) . iMakind 
( Il997h . and llwasawa et alj ( l2006h . Much of this work on SMBHs 
was based on earli er studies of stellar -mass black holes in glob- 
ular c lusters. Sig urdsson & HernquisJ ( 11993 ) and Kulkami et al.l 
( Il993r) considered the evolution of ~ 100 stellar mass black holes 
in globular clusters. They concluded that after mass segregation, 
most of these black holes are ejected out on a short time scale, 
and the globular cluster is left with none or a few black holes. 
Mass segregation and associated effects of stellar-mass black holes 
in a galactic nucleus with a central SMBH was a lso considered 
JMiralda-Escude & Goul j|2000l : lFreitag et alj2006h . 

The possible formation of systems with multiple SMBHs due 
to successive galactic mergers arises naturally in any model de- 
scribing the hierarchical assembly of galaxies. One approach to 
modeling SMBH growth involves constructing semi-analytic pre- 
scriptions of various characteristic processes, like mergers of galax- 
ies, formation of spheroids, star formation, and gas thermodynam- 
ics, coupled with merger trees of dark matter haloes. This ap- 
S iroach has been adopted, for example, bv lKauffmann & Haehneltl 
20001) ■ who also extended it to study possible formation of 
multiple SMBH systems and implications for the M(,h-o" rela- 
tion and density profiles obse rved in luminous ell iptical galaxies 
(Haehnel t & Kauff^maniill2002h . Another study bv IVolonteri et al.l 
(i2003i) followed merger trees of dark matter haloes and their com- 
ponent SMBHs using Monte Carlo realizations of hierarchical 
structure formation in the ACDM cosmology. They modeled dark 
matter haloes as singular isothermal spheres and calculated the in- 
spiral of less massive halos in more massive ones by using the 
Chandrasekhar formula for dynamical friction. Gas accretion to 
the SMBHs was modeled so as to reproduce the empirical Mbh~ 
cr relation and the SMBH dynamics was described with analytic 
prescriptions. In particular, the coalescence time of hard SMBH 
binaries was calculated from a set of coupled differential equa- 
tions based on scattering experiments involving the ejection of stel- 
lar mass from the loss cone due to the hard SMBH binary and 
the r esultant change in the hardening rate JOuinlanll 19961 : iMerritj 
I2OOCT) . For galaxies that underwent another major merger before 
their constituent binary SMBH coalesced, a three-body interac- 
tion was implemented between the binary and the intruder SMBH. 
They found that the smallest SMBH was kicked out of the galaxy 
in 99% of cases, while the binary escapes the galaxy in 8 % of 
cases. Thus, a significant fraction of galactic nuclei could end 
up with no SMBHs or offset SMBHs with mass lower than that 
expected from the Mbh-cr relation. These results were later ex- 
tended to incorporate recoil in the SMBH merger remnant due to 
asymmetric emission of gravitational waves, which mainly affected 
the Mbh-g' rela tion for low-mass haloes by increasing the scatter 
jVolonteri & Re es 2006; Volonteri 2007: Blechaet al.ir201lh . Sim- 
ilar semi-analytic models were studied by several other authors 
to understand the assembly of z ~ 6 quasars. However, most of 
these models ignored the dynamics of multiple SMBHs and as- 
sumed prompt coalescence fHaiman & Loeb 1999: Wvith e & Loebl 
2003b: Yoo & Miralda-Escude 2004: Tanaka & HaimarJ l2009l : 
Shen||2009). As a result, these models did not treat systems with 



multiple SMBHs. 

Lastly, SMBH assembly has also been studied using smooth 
particle hydrodynamic simulations that attempted to calculate ef- 
fects of both the gas physics as well as the gravitation al dynamics of 
the large-scale structure within and around galaxies faopkins et al.l 



ever, due to poor mass resolution and particle smoothening, these 
simulations cannot accurately calculate the detailed dynamics of a 
multiple SMBH systems. Indeed, in most of these studies, black 
hole coalescence occurs on scales smaller than the smoothening 
length, which is usually much larger than the expected separation 
of a hard SMBH binary. As a result, SMBH coalescence is imple- 
mented via a subgrid model. Here, we explore for the first time 
numerical simulations that incorporate the cosmological process of 
galaxy mergers in the cosmological context along with an accurate 
treatment of black hole dynamics. 



3 FORMATION OF MULTIPLE-SMBH SYSTEMS 

Unless they coalesce rapidly, or get kicked out of the host galactic 
nucleus, we expect multi-SMBH systems to form in galactic nuclei 
at high redshift due to mergers of galaxies if the typical black hole 
coalescence timescale is longer than the feeding timescale of new 
incoming black holes. In this section, we establish a simple theoret- 
ical framework for this formation path using analytical estimates of 
its relevant timescales: (i) the major merger time scale of galaxies: 
(ii) the time scale on which a satellite galaxy sinks to the center of a 
host galaxy so that a close interaction between SMBHs can occur: 
and (Hi) the time scale of SMBH coalescence. 



3.1 Time scale of incoming SMBHs 

iFakhouri et al.l ( 120101) have quantified the average merger rate of 
dark matter haloes per halo per unit redshift per unit mass ratio for 
a wide range of halo mass, progenitor mass ratios and redshift. The 
resul t is given by a fi tting formula derived from the Millennium 
( Spr ingel et alj 120050 and Millennium-II jBovlan-Kolchin et al.l 
120091) simulations: 
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d^dz 
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^exp 
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(3) 



Here, M is the halo mass at redshift z, and ^ is the mass ratio of pro- 
genitors. Mergers with ^ > 0.3 are considered major mergers. The 
best fit values of various parameters are a = 0.133, /? = -1.995, 
7 = 0.263, 77 = 0.0993, A = 0.0104 and | = 9.72 x 10^1 The 
average major merger rate per unit time is then given by 



dPL 
dt 



(M,z) 



I 



0.3 d^dz dt 



(4) 



IFakhouri et al.l J201(t) also provide a fitting formula for average 
mass growth rate of halos that can be used to calculate the halo 
mass at redshift z for use in equation ([3j, 



M(z) 
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(5) 



Using equation ^ we can now define the time scale of major merg- 
ers for a halo as 



'mii' 



dN,„ 



dt 



(6) 



l2006l : ISiiacki et al.l2007l : lLi et alj2007l : lHopkins et al.l2007l) . How 



The behavior of this quantity is shown in Figure [T] for three halo 
masses that discussed here: a Milky Way-like halo that has a mass 
Mo = 10'^ Mq at z = 0, the typical halo today that has mass Mq = 
lO"* Mo at z = 0, and rare haloes with mass Mq = 10"' Mq at 
z = 0. (In this paper, Mq always denotes the halo mass at redshift 
z = 0. We also refer to the average mass of such haloes at other 
redshifts, by e. g. M-4 and M-e- A halo with Mo = 10'^ Mq will 



© 201 1 RAS, MNRAS 000, 000-000 



4 Kulkarni & Loeb 



100 



10 






0.1 



0.01 



: 




1. 


7^^^ ^ > 0.3 


T- — 


19 ^ ^\ ^\ - 


Mo-10 Mq .^\ . . 
Mo=lolX '^X ^ 


" 


1 1 



0.1 



10 



Figure 1. Halo major merger time scale (mass ratio > 0.3), according to 
equation ^6), for haloes with mass Mq = lO'^ Mq (blue solid line), lO'** 
Mo (blue dashed line) and lO'^ Mq (blue dot-dashed line). The Hubble 
time is shown by the solid red curve. Major mergers are more frequent at 
higher redshifts. On average. Milky Way-sized haloes are not expected to 
undergo a major merger for z < 1 . Galaxy major merger time scale is always 
longer than the subsequent dynamical friction time scale. 



have M-_ 



2 X 10'" Mq. a halo with Mq 



10'^ Mo will have 



Mj=6 = 5 X 10" Mo. A cluster-size halo, with mass Mq = 10'^ 
Mq will have M~(, = 5 x 10'^ Mo and is expected to hold a single 
galaxy at that redshift.) This is the time scale at which we expect 
new (satellite) haloes to enter the halo. As expected, halo mergers 
are more frequent at higher redshift. At redshift z < 1 the major 
merger time scale for a Milky Way-like halo is greater than the 
Hubble time. 

After two dark matter haloes have merged, the smaller halo 
becomes a satellite halo within the virial radius of the host halo. 
It then takes this satellite a dynamical friction time to sink to the 
center of the host halo, so that the constituent galaxies can merge. 
As a result, the timescale for major mergers of galaxies is expected 
to be different that the time scale for major mergers of dark matter 
haloes calculated in Equation ^. 

The dynamical friction time scale is often estimated us- 
ing Chandrasekhar's formula JChandrasekhaill 19431 : iLacev & Cold 
ll993l : lBinnev & Tremainell2008h : 



fdf 



/dfQorb Mh. 



^dyi 



(7) 



In A Msa, 

where Mho,, and M^m are the masses for the host and satellite 
haloes respectively. In A is the coulomb logarithm, ©o^b is a func- 
tion of the orbital energy and angular momentum of the satellite, 
/df is an adjustable parameter of order unity and fjyn is the halo 
dynamical time scale calculated at the virial radius. Equation 
is valid only in the limit of small satellite mass in an infinite, 
isotropic and homogeneous coUisionless medium. Still, it has been 
used in the literature even for large satellite masses by modifying 
the Coulomb logarithm. In recent years, deviations from predic- 
tions by equation Q have been reporte d in both the M^a, <K Mt^gs t 
and Msa, <, Mhpst regim es (Taffoni e t al|2003l: [Monaco et aJjIlOOTl : 
iBovIan-KoIchin et al.''2008l; IJiang et alj2008| - IWetzel et al .ll2009l) . 
To correct the problems associated with Chandrasekhar's for- 
mula, several groups have developed full dynamical models of 
evolution of merging haloes JTavlor & BabulluOOll : lGnedinll2003l : 



iTafi^oni et al1l2003l ; Izeiitner et al.ll2005l) . For example, one of the 
approaches to overcome the limits o f Chandrasekhar' s formula is 
the theory of linear response (TLR; IColpi et al.lll999f) . TLR cap- 
tures the backreaction of the stellar distribution to the intruding 
satellite by correlating the instantaneous drag force on it with the 
drag force at an earlier time via the fluctuation-dissipation theorem. 
Tidal stripping of a satellite halo is an important ingredient in this 
formulation. In a singular isothermal sphere with ID velocity dis- 
persion (T and density profile p(r) = (r-/[2nGr^], TLR predicts a 
sinking time 



fdf =1.17 



GM,,t\nA 



(8) 



where e is the circularity (defined as the ratio between the angular 
momentum of the current orbit relative to that of a circular orbit of 
equal energy), rd, and Vcii are the initial radius and velocity of the 
circular orbit with the same energy of the actual orbit, and Ms is the 
mass of the incoming satellite halo. Nu merical simulations suggest 
a value of 0.4 - 0.5 for the exponent a Jvan den Bosch et al.ll 19991 : 
IColpi et alji999l ; IVolonteri et al.ll2003l) . 

Given the limitations of analytical treatments, we turn to re- 
sults of numerical simulations to understan d the dynamical fric- 
tion time scale. Using N-body simulations. iBovlan-Kolchin et al.l 
(2008) give a fitting formula that accurately predicts the time-scale 
for an extended satellite to sink from the virial radius of a host halo 
down to the halo's centre for a wide range of mass ratios and orbits 
(including a central bulge in each galaxy changes the merging time 
scale by < 10 %). Their fitting formula is given by 



'dyn 



r 



ln(l + 1/f) 



exp 



jcn(E) 



r^„(E) 
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where A = 0.216, b = 1.3, c = 1.9 and d = 1.0. Here f is the mass 
ratio Msat/Miiost, j is the specific angular momentum of the satellite 
halo, and jd, is the specific angular momentum of a circular orbit 
with the same energy E. This formula is expected to be valid for 
0.025 < ^ < 1.0, and for circularities Q3 < r) = jlj^AE) < 1.0. 
Mo st likely value of circularity in dark matter simulations is 77 a 
0.5 ( lBensoij|200^ ; IZentner et"S]l2005l:lKhochfar & Burkeril2006h . 
Lastly, it is valid for range of orbital energy -0.65 < raAE)lr^\^ < 
1.0. This covers the peak value of distribution seen in cosmological 
N-body simulations. We fix r^.iy(E)/r^ir = 1.0 and 77 = 0.5, which 
are the typical values found in simulations. 

We can now obtain the instantaneous merger rate of galax- 
ies by combining the halo merger rate and dynamical friction time 
scale. We closely follow the method of ^Sheg d 200ft) and write 



B,dM,^,z) = B[M,^,Ze(z,0] 



dZe 

Iz' 



(10) 



where B(M, ^, z) (per unit volume per unit mass per unit redshift 
per unit mass ratio) is the instantaneous merger rate of halos with 
mass M, progenitors with mass ratio ^ at redshift z, B„^\ is the same 
quantity for galaxies. The redshift Zeiz,^) is a function of z and ^, 
and is given implicitly by 

t{z)-t{Ze) = t^.M'Ze\ (11) 

where f(z) is the cosmic time at redshift z- IShenI ( 120090 finds that 
dzeldz is almost constant at all redshifts for ^ = 0.1 - 1 and can be 
approximated by 

dZe 



dz 



l-t-0.09[^'Mn(l + I/f)]- 



(12) 



for the fitting formula in equation ([9}. We assume this form in our 



© 201 1 RAS, MNRAS 000, 000-000 



Multiple Black Holes 



100 



Sh 


10 


p^ 




O 




m 




Q) 




a 


1 


o 




tn 




0) 




fi 




H 


0.1 




0.01 



Figure 2. A comparison between the feeding time scale of incoming black 
holes fi„ (black solid line; Eq. 1131 and the time scale of black hole coales- 
cence ?coal (black dot-dashed line; Eq. l23t . for a halo mass Mq = lO'^ Mq 
and considering only mergers with a mass ratio f = 0.4. The coalescence 
time t^a-ji has only a weak dependence on redshift because its dependence 
on Mbh and o" cancel out due to the Mbh-o" relation. This figure shows that 
at high redshift new black holes would arrive to the center of a galaxy faster 
than they could merge via dynamical processes. 



calculations. Once we have calculated Sgai(M, ^, z), we normalize it 
by n{M, z), the abundance of halo es of mass M at redshift z- We use 
the Sheth-Tormen mass function dSheth & Tormenll 19990 to calcu- 
late n(M, z)- This gives us the galaxy merger rate per halo per unit ^ 
per unit redshift, which is the galaxy's counterpart of equation Jsj, 
and which we denote by cWg^i/rfz. The rate of mergers of galaxies 
is the rate at which new black holes are added to the host halo's 
nucleus. Thus, the time scale of incoming black holes is 



fin 



rfA'gai dz 



dz dt 



(13) 



The result is shown by the solid black line in Figure [2] for a mass 
ratio of f = 0.4 and a halo that has mass of 10'' Mq at z = 0. 



3.2 Binary SMBH coalescence time scale 

In order to find whether there is a generic possibility of formation 
of systems with multiple SMBHs, we compare the time scale on 
which new black holes are added to the galactic nucleus at a certain 
redshift with the coalescence time scale of a binary SMBH at that 
redshift. 

As described in 21 ths formation and coalescence of a black 
hole binary is expected to take place in three stages. We define the 
coalescence time as the time that the binary spends in the second 
of these stages, that is the time from when the binary separation 
is a = fl/,, defined in equation l[T), up to when the separation is 
a = flgr at which point the binary enters the third stage of evolution, 
and gravitational waves become the dominant mechanism of en- 
ergy loss. For a hard binary, the dominant channel through which 
energy is lost is three-body interactions in which stars passing in 
close proximity to the binary are ejected at a much higher velocity 
Vcj = [GM,o,/a]''-, where M,o, is the total mass of the binary. The 
hardening time scale was quantified for a fixed stellar distribution 



bv lOuinlanI ( Il996l) . who found a time scale of 



ti,(a) ■: 



GpaH 



(14) 



where a is the binary separation, p and a are the density and one- 
dimensional velocity dispersion of the stellar background, and H 
is a dimensional parameter whose value was found from scatter- 
ing experiments to be 16 for a hard, equal-mass binary. In practice, 
however, the above expression for /;, is valid only during the initial 
stages of the binary's evolution. As the binary shrinks further, it 
ejects stellar mass from the central regions and modifies the stellar 
density p that appears in equation l ll4t . This fee dback can be qu an- 
tified using a simple analytical model given bv iMerritll 1 120001) , in 
which the binary evolution is described by two coupled equations, 
the first describing the binary's hardening due to the presence of 
stars. 



dt \a a 



(15) 



and the second describing the change in stellar density due to ejec- 
tion of mass by the hard SMBH binary, 



dMei 



dln(l/a) 



iM,„ 



(16) 



where Mej is the ejected mass , and J is anot her dimensionless pa- 
rameter that was measured bv lOuinlanI ( 11996) to be close to unity 
and nearly independent of a. 

By assuming a singular isothermal sphere profile for the stel- 
lar density and assuming that the ejected stellar mass c auses a 
consta nt-density core to form at the center of this profile, IMerrittl 
( 120000 finds that evolution of the binary separation can be described 
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where a;, is as defined in Equation ([T}, a(finii) = a/,, and ?o is given 
by 

97ry^ /Mie 
H \ 2^2 / ^ 
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(18) 



This result is found to closely match with the evolution observed in 
N-body simulation. 

On the other hand, the timescale for emission of gravitational 
waves is given by 



(19) 



256 G^mimiM,^ 



As a result, the binary will continue to harden only up to the time 
when hardening time f;, = fgr, after which it will coalesce rapidly 
due to gravitational wave emission. Using equati on ( I17t . it can be 



equatK 
( iMerri 



^'-AllnAr, 



shown that this occurs when a = Ogr where ( lMerritli200Qi) . 
and 

xO.2 / ,, \0.4 



m2 1 \ 21112 I \ t / 



(20) 



(21) 



Here mj and m2 are masses of the components of the SMBH binary. 
Finally, we can again use equation dlVt to calcu late the time i t takes 
for the binary to shrink from a = «/, to a = a^ ( IMerri tJbOOOf) : 



fcoai~8foA"'|lnAr 



8/5 



(22) 
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which can be simphfied as 



f,„al«1.4xl0'« 



/m2\°'^/Mtot 
[mil \2m2 



M„ 



10«Mo/\200km/s 



(23) 



Clearly, there is a possibility for the formation of multiple- 
SMBH system if f;,, < fcoai- These two time scales are compared 
in Figure|2]for a halo that has a mass of Mo = 10'^ M© at z = 0. 
For simplicity, we have fixed the mass ratio of merging haloes to 
be ^ = 0.4. At each redshift, we calculate fi„ from equation l ll3b . 
In order to estimate fcoai at a given redshift using equation l |23l ), 
we first infer the mass of the halo at that redshift from the fitting 
function for the halo's assembly history from equation (O. We then 
assume that a galaxy belonging to a satellite halo with mass ratio ^ 
has merged with this host halo at this redshift. 

In order to estimate the mass of bla ck holes in the nuclei of 
these galaxies, we follow the approach of IHofi^man & Loebl 1 12007!) 
in employing the Mbh-o" relation. The virial velocity (defined as the 
circular velocity at virial radius) for a halo of mass M at redshift z 
is given by 



. = 23.4 



M 



lO^h-'Me 



1/3 



n,„ A, 



Q.;„ 18;r2 



1+z 
10 



1/2 



km/s. 



where 
Q.- = 



n,„(i + z)3 



(24) 



(25) 



a,(l+z)3+fiA + fit(l+z)^' 
and Ac is the overdensity of the halo relative to the critical density, 
given for the ACDM cosmology by 

A, = 18;r2 + 82d-39d^ (26) 

where d = £!;„ - I JBarkana & LoebluOOlh . Further, we equate the 
halo virial velocity with the circular velocity vv of its constituent 
spheroid an d obtain the vel ocity dispersion of the spheroid using 
the relation ( lFerraresel2002l) 



314 



km/s. 



208km/s 
This combined with the Mbh-o" relation (ITremaine et al.ll2002h 

208km/s ~ 1.56 x IO^Mq 



(27) 



(28) 



gives 



Ml 



halo 



, lO'^Mp 



8.28 



Mhi 



10** Mf 



Q,„ A, 



n-l/2 



n~, I8;r2 



{l+z)-"\ 



(29) 



We obtain the black hole masses in the host and the satellite 
haloes using equation ( |29l > and use the spheroid velocity dispersion 
from equation J27| ( to estimate the coalescence time from equation 
([23}. The result is shown by the dashed line in Figure|2] 

At high redshift, early on in the assembly history of a halo, 
the galaxy merger rate is higher than the SMBH binary coalescence 
rate and systems with multiple SMBHs can form. Note that the time 
scale fcoai obtained above will change if effect of loss-cone_replen- 
ishment and gas are taken into account. However. I Yul ( 120021) finds 
that in realistic spheroidal galaxies, even loss-cone replenishment 
is insufficient to cause early coalescence. 



4 EVOLUTION OF MULTIPLE SMBHS 

We have described the literature on systems with more than two 
SMBHs in ^ If the infalling SMBH is less masssive than either 




Figure 3. An example merger tree form the Millennium .simulation of a 
halo that has a mass of Mq ~ lO'^ Mq. This plot shows major mergers 
(mass ratio > 0.1) in all branches of the halo's merger tree. 



of the components of a pre-existing binary then we expect the ulti- 
mate o utcome to be eject i on of the smaller SMBH and recoil of the 
binarv. lHoff^man & Loebl ( 120070 studied the statistics of close triple 
SMBH encounters in galactic nuclei by computing a series of three- 
body orbits with physically motivated initial conditions appropriate 
for giant elliptical galaxies. Their simulations included a smooth 
background potential consisting of a stellar bulge and a dark matter 
halo, and also accounted for the effect of dynamical friction due to 
stars and dark matter. They found that in most cases the intruder 
helped the binary SMBH to coalesce via the Kozai-Lidov mecha- 
nism and by scattering stars into the binary's loss cone. In this case, 
the intruder itself was left wandering in the galactic halo, or even 
kicked out of the galaxy altogether. It was also found that escape of 
all three black holes is exceedingly rare. 

Dynamical evolution of multiple massi ve black holes in glob- 
ular clusters has received mu ch attention dKulkarru et al.l 1 1 9931 : 
ISigurdsson & Hemquistlll993h . From these studies, it is expected 
that systems with more than two SMBHs will last for about a cross- 
ing time. 



5 SIMULATIONS 

In order to accurately calculate the formation and evolution of 
galactic nuclei with multiple black holes, we perform direct- 
summation N-body simulations of galactic nuclei merging in a cos- 
mological context. This essentially involves generating physically 
consistent initial conditions for galactic nuclei with SMBHs at high 
redshift and evolving them while taking into account the mergers 
of such nuclei and the resultant close interaction of their SMBHs. 

We obtain merger histories of galactic nuclei by extracting 
merger trees of gravitationally bound subhaloes from the Millen- 
nium Simulation Databasa4 which stores results of the Millennium 



Simulation dSpringel et al.ll2005f) . The Millennium Simulation is a 
pure dark matter simulation with a ACDM model with 2160' par- 
ticles in a periodic cube 500 h^'Mpc on a side. This corresponds 

^ http://www.mpa-garching.mpg.de/millennium/ 
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Simulation 


Mass of halo at z = 


(Mo) 


Max. BH no. 


SMBH Coalescences 


SMBH Escapes 


LI 


1.21x10"* 




4 


7 


2 


L2 


1.31 xlO''* 




2 


1 


1 


L3 


1.31x10''* 




2 


3 


2 


L4 


1.24 xlO'"* 




2 


5 


5 


L5 


1.28 X 10'"* 




5 


8 


4 


L6 


1.31 X 10"* 




6 


6 





L7 


1.23 X 10'" 




3 


2 





L8 


1.31 X lO'* 




2 


3 


1 



Table 1. Summary of simulations and results for haloes that have a mass of Mq ~ 10 Mq. The maximum BH number denotes the number of black holes in 
the biggest BH group found in a simulation. The last two columns show number of BH coalescences and escapes in the simulation. A halo with Mq = lO'* 
M© has average mass M,=6 = 5 X lO" Mq. 



Simulation 


Mass of halo at z = (Mq) 


Max. BH no. 


SMBH Coalescences 


SMBH Escapes 


HI 


1.25x10'^ 


6 


4 


3 


H2 


1.65 X 10'^ 


2 


1 


1 


H3 


1.81xl0'5 


3 


2 





H4 


1.24xl0'5 


5 


6 


3 


H5 


1.37 X 10'^ 


3 


7 


1 


H6 


1.40 X 10'^ 


4 


3 





H7 


1.41 X 10'5 


6 


9 


1 


H8 


1.45xl0'5 


3 


4 


1 


H9 


1.46 X 10'5 


2 


2 





HIO 


1.48 X 10'^ 


4 


7 


1 


Hll 


1.54 X 10'^ 


2 


1 


1 


H12 


1.59xl0'5 


5 


10 


1 


H13 


1.66xl0'5 


8 


15 


4 


H14 


1.71x10'^ 


4 


3 





H15 


1.81 X 10'5 


4 


20 


7 


H16 


1.86x10'^ 


3 


7 


4 


H17 


4.04 X 10" 


8 


11 


2 



Table 2. Summary of simulation runs with haloes that have mass Mq > 10 M© at z = 0. Various columns are same as Table[T] A halo with Mq = 10 M© 
has average mass M~^ = 5 X lO'^ M©. 
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Figure 4. Evolution of single and binary SMBHsin our simulations, (a) The left hand panel shows evolution of the .v-component of the position of a 9.95 X 10 
M© back hole near the centre of a Hemquist bulge of mass 5.41 X 10^ M© and scale length of 0.2 kpc. The particle mass is 5.41 1 X 10^ M©. The secular motion 
is due to that of the cusp, (b) The right hand panel shows evolution of the separation between SMBHs in a binary with initial separation 2 kpc and eccentricity 
0.5. The black hole masses were 8.65 X 10 M© and the binary evolved near the center of a Hernquist halo with mass 5.41 X 10 M© and scale length of 10.0 
kpc. The particle mass is 5.41 1 X 10^ M©. 
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to a particle mass of 8.6 x 10** h"' Mq. The output of this simula- 
tion is stored in 64 snapshots between z = 127 and z = 0. Particles 
in each snapshot are grouped into friends-of-friends (FOF) clus- 
ters that are expected to correspond to virialised structures. Each 
FOF halo contains substructure of gravitationally bound subhaloes 
that can be related to each other across snapshots as progenitors 
and descendants. Because a halo can contain multiple galaxies, we 
expect the subhalo merger tree to reflect the merger history of the 
galaxies within a halo. Since the goal of this paper is to under- 
stand formation and evolution of systems of multiple black holes 
due to the hierarchical merger history of a galaxy, we extract sub- 
halo merger trees from the Millennium Simulation Database. Each 
such merger tree typically shows growth of a subhalo via accretion 
of dark matter particles and via mergers. We process these merger 
trees to keep only major mergers, which we define to be mergers 
having mass ratio larger than 0. 1. To identify the mass ratio of two 
subhaloes, we use the masses of the distinct FOF haloes that these 
subhaloes were a part of before the FOF haloes merged. This is to 
account for the mass loss of the satellite subhalo due to tidal strip- 
ping after it enters the FOF group of the host subhalo, but before 
the eventual merge r of the two subhaloes. (See discussion in §5 of 
iBundv et alj|2007l .) Figure |5] shows the resultant merger history of 
a Milky Way sized halo. The main reason behind removing minor 
mergers from our calculation is that for such mergers the dynami- 
cal friction time taken by the satellite halo to reach the center of the 
host halo is longer than the Hubble time. As a result, in such merg- 
ers, we do not expect the constituent galactic nuclei of these haloes 
to interact closely. Since, as we describe below, we model only the 
spheroidal galactic nuclei in our simulations, we only need to ac- 
count for mergers in which such nuclei will closely inte ract. This 
approach is very similar to that used by iLi et alj ( 120070 ■ with the 
main difference being our use of direct-summation N-body simula- 
tions instead of SPH simulations. 

Once we have a galaxy merger tree, we set up the initial condi- 
tions of our simulation in the "leaves" of the tree, that is, in haloes 
that do not have a progenitor, and follow the evolution using an N- 
body calculation. The initial conditions of our simulation consist of 
a stellar spheroid with a Hernquist density profile. 



p(r): 



M 



2n r(r + aY 



(30) 



where M is the total mass of the spheroid and the scale length a is 
related to the half mass radius tj/t of the spheroid by a = 0.414ri/2. 
Values for the p arameters M and a we re obtained from the halo 
mass as follows JHoflFman & Loebll2007n . We first obtain the black 
hole mass M(,h from the halo mass M^aio using Equation J29I ). 
We then use the empirical relation between the SMBH mass and 
the spheroid's virial mass JMagorrian et alJll998l : lMarconi & Hung 



the spheroid s virial m ass ( IMagornan et alJ lr 
l2003l : |Pengetal.ll2006h to obtain the latter as 



M, 



,ph 



: 4.06 X lO'^M, 



Mh, 



lO'^Me 



(31) 



The virial mass of the spheroid is related to its velocity dispersion 
(Tj, and half light radius R^ by 



M,pi 



kR,o-l 



(32) 



We follow iMarconi & Hunt! ( 120031) and set fc = 3 to get an aver- 
age ratio of unity between this mass estimate and the dynamically 
measured masses of galaxies. The velocity dispersion in the above 
equation is usually measured over either a circular aperture of ra- 
dius /?j/8 or a linear aperture of length R^. These two methods are 



in essential agreement, as argued by iTremaine et alj ( l2002h . As- 
suming a constant mass-to-light ratio for the Hernquist profile, we 
ha\eRg = 1.815a and the velocity dispersion at radius ^^/S is given 
by 



0.104GM 



(33) 



This lets us obtain the value of the parameter M of the Hernquist 
profile as M = 1.765Msph. The scale length a is readily obtained as 



3/^1 <' 



(34) 



where o-f,i, is obtained using the M - cr relation of equation ( I28t . 
Having obtained a density profile for the bulge, we place a black 
hole at its center and set the black hole mass to be ten times that 
obtained from equation ( |29t . This factor of ten is introduced to 
keep the ratio between the black hole mass and the particle mass 
high enough dMilosavlievic & M erritt 200 ll: iMakino & Ebisuzakil 
Il996) . We confirm that the radius of influence r-i^f = Gm^^lcr^ of 
this black hole is still much smaller than the a. Velocities of the stars 
in the spheroid are then generated from the unique, isotropic ve- 
locity distribution that corresponds to the gravitati onal potential of 
the density profile in Equation ([30]l and the SMBH ( ITremaine et al.l 
12002) . These initial conditions are then scaled to standard N-body 
units of G = 1, M = 1 and E = -0.25, where M is the total mass 
of the system and E is its total energy (^Heggie & Mathieulll986l : 
lAarsethl2003h . In these units, in virial equilibrium, the mean square 
velocity (v~> = 1/2 and the system's crossing time is f^r = 2 ^/2, in- 
dependent of the number of particles. The conversion factors from 
physical units to these N-body units can be easily obtained via di- 
mensional analysis. 

Note that we ignore presence of gas in this set-up. Simu- 
lations of binary BHs in gaseous environment have not reached 
sufficient resolution to establish th e role played by gas in evolu- 
tion of SMBHs in g alactic nuclei JMerritt & Milosavlievi3l2005l : 
IColpi & Dottill20 09 ). Moreover, we expect that at high redshifts, 
quasar activity triggered by galaxy mergers could efficiently drive 
gas away from the shallow potential wells of the galaxies. 

To perform the actual dynamical evolution of this system, 
we use the direct- summation code NBODY6 written by Sverre 
Aarseth ( Aarsethl l 19991 . 1200 3). This code has been well-tested for 
various applications since around 1992. Its purpose is to perform 
an exact integration, without particle softening, of a large num- 
ber of particles. It integrates equations of motion of individual 
partic les using a fourth-order Hermite method with block time 
steps ( IMakino & AarsethI 19920 . T his integrator is coupled with the 
Ahmed-Cohen neighbour scheme lAhmad & Cohenll 19730 . which 
selects a subset of neighbours of a particle whose forces on it are 
calculated at a higher time resolution that other, more distant, parti- 
cles. This scheme reduces the computational cost from O(N^) to 
about (9(A''^). Close two-body encounters are treated using the 
Kustaanheimo-Stiefel (KS) regularization method that eliminates 
the r = singularity in Newtonian gravity by using a coordinate 
transformation. Triples, quadruples and compact subsystems of up 
to six particles (ca lled "chains") are treated using the chain regu- 
larization method jMikkola & Aarsethll 19901) . Details of the vari- 
ous algorithms in this code and their implementation are given by 
lAarsethl ( 120030 . In all simulations reported in this paper, the time- 
step parameter for irregular force polynomial, rjj, and the time-step 
parameter for regular force polynomial, rjf, are set to 0.02. The en- 
ergy tolerance is set to Qe = 4x 10"^ and the regularized time-step 
parameter is set to rju = 0.2. 
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We check the stability of our initial conditions by evolving 
standalone realizations of the Hemquist bulge with a central BH 
and then traverse the merger tree of a given halo using NBODY6, 
starting from the initial conditions as described above. We scale the 
physical time between two successive nodes of the tree to N-body 
units and run NBODY6 for that duration. If a merger happens at 
a certain node, we place the two galactic nuclei at a distance of 2 
kpc apart and evolve in an head-on approach. Although such head- 
on mergers would be unlikely, we choose it to reduce the compu- 
tational time while still retaining some realism. When two galax- 
ies, that are in equilibrium separately, merge we expect some tran- 
sient response in the resulting d ynamics. However, as discussed by 
iMilosavlievic & MerrittI ( 1200 ih . any such effects in the dynamics 
of the central regions of the merger remnant of these galaxies are 
essentially negligible. 

Under these conditions, the component black holes approach 
after a merger event and the remnant galactic nucleus is left with 
two black holes, which gradually harden due to dynamical friction 
and three-body interactions with stars in their vicinity. Black hole 
coalescence is implemented in our simulation by monitoring the 
separation of hard black hole binaries. Once members of a SMBH 
binary get closer than a fixed distance dau, we replace them with 
a single black hole with mass equal to the sum of the masses of 
component black holes. In all the runs reported in this paper, we set 
c4,it = 0.1 pc. Note that this is the only mechanism in which black 
holes grow in our simulations. Thus, the initial SMBH masses are 
set according to the M - a relation, but the later growth of these 
SMBHs occurs only via coalescence. 

Recoil due to anisotropic emission of gravitational waves is 
a natural consequence of asymmetric merger of bla ck holes, ei- 
ther due to unequ al masses or due to unequal spins jPerejl 19621 : 
iBekensteinl 1 1 9731) . Until recently, it was unclear whether this re- 
coil is large enough to be astrophysicaly relevant. However, recent 
results from numerical relativity have revealed t he resultant kick 
velocities in a va riety of merger configurations jPretoriusI 120051 : 
iBaker et alj2006l) . When the black hole spins are aligned with each 
other and with the orbital spin, these simulations find revoil ve- 
locity o£j;recoU_^_200^ms_^_jBaker_e^al. 2006; Gonzalez et al.l 
l2007l : lHerrmann et alj2"007l : lLousto & Zlochower.200^ . In the ab- 
sence of spins, this recoil velocity is only a function of the ratio of 
black hole masses. For random orientations of s pins, recoil veloc- 
ities as hi gh as 2000 km s~' have b een obtained JCampanelli et al.l 
l2007al lbl). iBogdanovic et alj ( 120071) argue that a circumbinary gas 
disk can align the binary spins with the orbital axis thereby reduc- 
ing Vrecoii to about 200 km s"' . In our simulations we assume a con- 
stant kick velocity of 200 km s"', which we impart to the remnant 
of every unequal-mass binary SMBH coalescence. 

We follow the approach of lMakino & Aarsethl ( 119920 and keep 
the particle number fixed at A^ = lO'* throughout the simula- 
tion. Thus, at every merger, we combine particles in each merg- 
ing galactic nucleus and double the particle mass. This lets us 
keep the particle number high throughout the merger tree of the 
halo. The ratio of black hole mass to the stellar mass is typically 
a few hundred, which is also roughly the ratio of the spheroid's 
total mass to the black hole's mass. These values are compara- 
ble to other simulations of th is kind JMakino & Ebisuzakj|l996l : 
IMilosavlievic & Merrittll200ll) . 

In summary, the unique features of our simulations are: (i) 
kinematically consistent initial conditions with black holes; (ii) cal- 
culation of mergers of galactic nuclei in a cosmological setting us- 
ing merger trees extracted from cosmological N-body simulations: 
(in) calculation of merger of galactic nuclei resulting in a formation 
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Figure 5. Number of black holes as a function of redshift in a simulation 
with Mo = 1.29 X lO'" Mq. 



of SMBH binaries starting from the results of each nucleus having 
evolved in isolation; and (iv) accurate calculation of SMBH-star 
and SMBH-SMBH dynamics throughout the assembly history of a 
galactic nucleus and its constituent SMBH with the effect of gravi- 
tational wave recoil taken into account. 



6 RESULTS 

We perform some basic checks on our code, such as ensuring en- 
ergy conservation and stable evolution of equilibrium systems. In 
all of our runs, the relative error in the total energy is maintained at 
\^E|E\ < 4 X m-\ The treatment of BH-BH and BH-stai' interac- 
tion is handled by the original nbody6 code, and is expected to be 
accurate. One caveat here is that the neighbour criterion in nbody6 
for regularization of close particles is based on inter-particle dis- 
tance. As a result, while evolving a set of particles in the vicinity of 
a massive BH, the code either selects a large number of particles for 
chain regularization, or selects every close pair of particles for two- 
body regularization. This usually slows down the code. Indeed, in 
three of our runs the code run time exceeded practical constraints 
because of this effect. These three runs are excluded from the re- 
sults presented below. 



6.1 Dynamics of single and binary SMBHs 

In a stellar environment, a single SMBH exhibits a random fluctuat- 
ing motion arising due to discrete interactions with individual stars. 
As a result, the effect of the stellar environment on the SMBH can 
be decomposed into two distinct components: (1) a smooth compo- 
nent arising due to the large scale distribution of the whole system, 
and (2) a stochas tic fluctuating p art coming form the interation with 
individual stars tChatteriee et ali 2002). This random motion is il- 
lustrated in the left hand panel of Figure[5] which shows evolution 
of the x-component of the position of a 9.95 x 10^ Mq back hole 
near the centre of a Hemquist bulge of mass 5.41x10^ Mq and scale 
length of 0.2 kpc. The particle mass is 5.41 1x10' Mq. As expected, 
the SMBH wanders around due to stochastic interactions with the 
stars in its vicinity. The mean square amplitude of these fluctuations 
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Figure 6. Histograms of ejection velocities of BHs. Left: Velocities of ejected black holes in all of our high mass runs. Note that this does not include ejected 
black holes with the highest velocities (> 2000 km s"'). Right: number of ejections as a function of redshift in our high mass runs. 



is exp ected to be JChatteriee et alj|2002l : iMilosavlievic & MerrittI 
l2003ah 



<x^> 



■2\.^. 



m, 2 

COI 

msH 



(35) 



where rcoic is the radius within which the stellar distribution flat- 
tens out. The Hernquist distribution that we have used here does 
not have a well-defined core, since the density keeps rising as r"' 
near the origin. Milosavlievic & Merritt (2003a) argue that the ef- 
fective core radius for such distribution can be taken as the ra- 
dius of influence of the black hole. The resultant mean square 
value of fluctuations is somewhat smaller that that for Figure |5]by 
roughly a factor of 2 as is known to happen in N-body simulations 
JOuinlan & Hemouist 1997 ; Milosavlievic & Merritt 200 3a). 

As described in fl the evolution of a binary black hole in 
a gas-poor galaxy takes place in three stages. Right hand panel 
of Figure [5] shows evolution of the separation between SMBHs in 
a binary with initial separation 2 kpc and eccentricity 0.5 in our 
code. The black hole masses were 8.65 x lO'' Mq and the binary 
evolved near the center of a Hernquist halo with mass 5.41 x 10' 
Mq and scale length of 10.0 kpc. The particle mass is 5.411 x 10^^ 
Mq. In the first stage of evolution, the SMBHs sink to the cen- 
tre of the galactic nucleus by losing energy via dynamical friction 
and become bound to each other. This stage ends when the sep- 
aration bet ween the SMBHs is equal to the radius of influence of 
the binary JMerritt & Milosavlievicl2005r) . In the second evolution- 
ary stage, the binary loses energy predominantly ejection of nearby 
stars via three-body interaction. The binary loses energy rapidly in 
this stage, which continues until t x 200 Myr for the case depicted 
in Figure |5] The final stage of the SMBH binary evolution begins 
when the rapid hardening of the second stage stops. This happens 
when the binary semi-major axis takes the value given by Equation 
0). The binary semi-major axis is related to the separation r by 

I 2 V- 

- = , (36) 

a r fi 

where v is the relative velocity of the BHs and u is the reduced mass 
('Makino & Funato" 2004 ISerczik et al.l 120061 ; iMerritt et al.ll2007l ; 
[Khan et al. 201 1). In A^-body simulations, the last stage is known to 
have a dependence on the number of part icles A^ such that the hard - 
ening rate decreases with increasing N dMakino & Funatdl2004h . 
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Figure 7. Escape velocities from the bulges of haloes in our three cate- 
gories of present-day masses of haloes. Solid line: Mq ~ lO'^ Mq, Dashed 
line: Mq ~ lO'"* Mo, Dot-dashed line: Mq > lO'^ Mq. Note that these are 
average values computed from the fitting functions to the Millennium sim- 
ulation. Therefore, case by case comparison with our runs is not straight- 
forward. 



For real spherical galaxies, the binary separation would stop evolv- 
ing after this point because the loss cone is empty. 



6.2 Evolution of nuclei with multiple SMBHs 

We now run the simulation along merger trees of haloes drawn from 
the Millennium simulation as described in Section[5] These simula- 
tions are described in Tables[T]and[2] We randomly select 8 haloes 
with mass Mg around 10'* Mq at z = 0. These correspond to the 
typical haloes (M « M.) in the present epoch. We also randomly 
select 17 haloes whose present-day mass Mq is in excess of 10''' 
Mq. These are rare, high mass haloes t hat are expected to host the 
redshift 6 SDSS quasars JLi et al.l2007r) . Additionally, we have also 
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Figure 8. Number of black holes as a function of redshift in a few of our simulation runs. 




Figure 9. Projected stellar density contours in the presence of a binary in the simulation H5. Each panel is 400 pc on a side. Clockwise from top left to bottom 
right, the redshifts are z = 10.073, 8.54, 7.27, 6.19, 5.28, and 4.52. The total time span is about 800 Myr. Core-SMBH oscillations are clearly visible. 
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Figure 10. Projected stellar density contours in the presence of multiple BHs in the simulation H4. Each panel is 100 pc on a side. The total time span, 
clockwise from top left to bottom right, is about 1 Gyr. Most BHs are stripped of their cusps in nuclei with multiple BHs. 
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Figure 11. The fraction of runs with multiple SMBHs at different redshift bins for haloes with a mass Mq ~ IO'^Mq at j = 0. The results of these runs are 
summarised in Table [T] The three panels from left to right describe the occurrence of systems with more than 2, 3 and 4 black holes respectively. At each 
redshift, this number can be interpreted as the likelihood of finding such systems in haloes of mass Mq ~ 10''' M© at z = 0. It is seen that systems with multiple 
SMBHs ai-e rare at redshift z < 2. Note that a halo with Mq = lO''' M© will have Mj=6 = 5 X lO" Mq. 



simulated 11 haloes with present-day mass similar to the Milky 
Way halo (Mq ~ 10'" MQ).Using the prescriptions described in the 
previous section, and using the N-body integrator, these simulations 
tell us about the effect of multiple mergers of galactic nuclei with 
SMBHs. 

Figure [5] shows results from a typical simulation run, for a 
halo of mass 1.29 x 10''' Mq. We plot here the number of BHs in 
the bulge in the main branch of the galaxy's merger tree at vari- 
ous redshifts. It is seen that the central bulge has more than one 
SMBH for a wide redshift range (2 < z < 6; about 2.5 Gyr). For 



3 < z < 5 (about 1 Gyr) the bulge holds more than 2 BHs. The 
maximum number of BHs interacting within the bulge in this sim- 
ulation is 6. Lastly, the number of BHs reduces to one well before 
z = due to coalescences and ejections. Note that at the highest 
redshifts (z > 6) there are no BHs in the central bulge. This is sim- 
ply an artifact of the limited numerical resolution of the Millennium 
simulation, because of which the halo merger tree is not resolved 
at these redshifts. To ensure that this does not affect our results for 
z < 6, we set up initial conditions at z ~ 6 such that the BHs are on 
the M-cr relation, and by using a Hernquist bulge with inner slope 
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Figure 12. The fraction of rans with multiple SMBHs at different redshift bins for halo masses Mq > lO'^M© at z = 0. The results of these runs are summarised 
in Table |2] The three panels from left to right describe the occurrence of systems with more than 2, 3 and 4 black holes respectively. At each redshift, this 
number can be interpreted as the likelihood of finding such systems in haloes of mass Mq > IO'^Mq at z = 0. It is seen that systems with multiple SMBHs 
are rare at redshift z < 2. These results can be compared with those in figure [TTI Nuclei with multiple SMBHs are more likely in high mass haloes because of 
higher merger rate. Note that a halo with Mq = lO''' Mq will have M-^(, ~ lO'^ Mq. 



-1. In the absence of gas, the systems with multiple SMBHs form 
generically, in high mass haloes with frequent of major mergers. It 
is evident than such systems are usually short-lived and most often 
these nuclei contain a single SMBH at z = 0. Most SMBHs es- 
cape into the halo, where they join a population of wandering black 
holes or escape the halo completely. 

Similar results from a few other simulation runs for haloes 
with mass Mq ~ 10'"* Mq at z = are shown in Figure[8] Most of 
these runs have features similar to the run described above. Multi- 
ple BH systems form generically and last for 2-3 Gyr. Importantly, 
most of these galaxies end up with a single SMBH in their central 
bulge. This is in contras t with expectation s from some simple ar- 
guments in earlier work ( IHut & Reeslll992') . About 5% of galaxies 
in our simulations end up with no BHs in their centres at z = 0. 
Tables[T]and|2]summarize these features of all our simulations. The 
last columns of these tables show the cumulative number of BHs 
that were ejected out of the galactic nucleus throughout the run ei- 
ther due to recoil associated with emission of gravitational waves 
or due to many-body interaction between the BHs. We find that 
for most triple and quadruple SMBH systems in our calculation, 
gravitational wave recoil is the dominant mechanism for SMBH 
escape. Many-body interaction between SMBHs was the dominant 
cause only when the number of black holes was more than four. 
Consequently, for low-mass galaxies in which the number of BHs 
is small, almost all escapes were because of gravitational wave re- 
coil. Whereas in our low mass galaxy simulations, larger number 
of coalescence usually results in large escapers, in the high mass 
galaxy simulations, coalescence often does not lead to escape. In 
high mass galaxies, BH-BH interaction is the dominant mechanism 
behind escaping SMBHs. Figure[6]summarizes this. The right hand 
panel shows that most ejections happen at high redshifts. Typical 
ejection velocities are seen in the left hand panel. Ejection veloci- 
ties are spread out up to 200 km s"' , which is the GW recoil kick in 
our simulations. Note that this plot does not show kicks with very 
high velocities, which we describe below. 

With the prescription that we have adopted in this paper, we 
find that SMBH coalescence happens in each one of our simu- 
lations. Tables [T] and |2] give the number of BH coalescences oc- 
curring in our simulations. Due to the limitation on the particle 



number, our simulations implement BH coalescence by replacing 
a bound binary BH by a single BH whose mass is equal to the to- 
tal mass of the binary. As an example. Figure |9] shows the merger 
of two bulges beginning from initial conditions at redshift 6.7 in 
the run H5. In Figure |9l the hardening radius is a/, = 0.5 pc at 
ti, = 500 Myr. We find the the BHs remain associated with their 
host cusps until cusp coalescence. It is known that by increasing 
the effective mass of the BHs, this increases the rate of coales- 
cence of the BHs by as much as ~ 6 times compared to the dy- 
namical friction time scale. We also see the homology of density 
structure before and after the merger, as reported previously in the 
literature d Milosavlievic & Merritll200ir) . However, one prominent 
difference from previous works is in the evolution of the density 
profile in the later stages of the merger. In our simulations, each co- 
alescence event is followed by recoil of the remnant at 200 km s"', 
which at high redshift, usually results in the escape of the SMBH 
from the galaxy. At relatively low redshifts, the recoiled SMBH re- 
turns to the nucleus in few hundreds Myr. Because of this recoil, 
the remnant BH is detached from its cusp immediately. At the re- 
coil speed implemented here, this happens at a much smaller time 
scale that the local crossing time scale. As a result, the only effect 
of the remnant on the cusp is due to subsequent core passages. 

Usually, most coalescences are assumed to take place due to 
BH hardening via BH-star encounters. In gas-free systems, this 
leads to the final parsec problem. In our simulations, we find that in 
high mass haloes, roughly half of the SMBH coalescences are due 
to three-body scattering with intruder SMBHs. This is expected, 
since in spite of higher major merger rate, high mass galaxies in our 
model are still left with at most two SMBHs at z = 0. The dominant 
mechanism of coalescence is then three body interactions. Figure 
[8] shows an example of the evolution of a multiple BH system that 
undergoes three coalescences due to BH-BH dynamics. We find 
violent oscillations of the cusp-BH system as shown in Figure |9] 
This has significant impact on the density distribution of the core, 
and also results in off-centre BHs, which slowly return to the centre 
of the cusp due to dynamical friction. 

About 10% of SMBH ejections in our simulations occur at 
very high speeds of > 2000 km s"'. In haloes with Mq ~ 10'^ 
Mq these SMBHs will linger in the outskirts of the halo for 2-10 
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Figure 13. Evolution of density profile for simulations H3 and H5 in N-body units. The solid line is the oiiginal Hernquist profile with an inner logarithmic 
slope of y « - 1 . Dashed line shows the profile after one SMBH binary coalescence, dot-dashed line after the second coalesence and the dotted line after the 
third coalescence. These plots are shown in N-body units to scale out the doubling of the half-mass radius. See text for details. 



Gyr as can be seen by comparing with the bulge escape speeds 
in Figure |7] The SMBHs in the wandering phase that are intro- 
duced via this mechanism have markedly different properties than 
the BHs introduced due to galaxies that have not yet reached the 
host g alaxy's center so as to have a close encounter i Volonteri et al.l 
l2003h . The main difference is that our ejected black holes are much 
more massive than those in the other category. Moreover, the veloc- 
ity of ejected SMBHs will typically be higher that black holes in the 
other category, which have already experience significant dynam- 
ical friction. Three of the 30 BH ejections in our runs are ejected 
binaries. 



6.3 Likelihood of nuclei with multiple SMBHs at high 
redshift 

From the results of our simulations, we can estimate the likelihood 
of galactic nuclei with multiple black holes at high redshifts. The 
histograms in Figures[TT]and[T2]show fraction of runs with multiple 
SMBHs at each redshift for haloes with present-day masses of ~ 
lO'"* Mq and ~ 10'^ Mq, respectively. The three panels from left 
to right describe the occurrence of systems with more than 2, 3 
and 4 black holes respectively. At each redshift, this number can be 
interpreted as the likelihood of occurrence of such systems at that 
redshift. 

Systems with more than 2 SMBHs are generically expected in 
the central galaxies of haloes with Mq > lO'"* M© at around z > 3. 
On the other hand, few galaxies hold multiple black holes at red- 
shifts z < 2 because the galaxy merger rate is low at these redshifts 
and the BHs have sufficient time to coalescence. This is consistent 
with the expectation from our heuristic analysis of Section [3] In 
other words, multiple black hole systems are numerous at around 
redshifts of 6, when there are many major mergers in the system. 
Our numerical simulations show that such systems can exist in suf- 
ficiently long-lived configurations of SMBHs separated on pc-kpc 
scale. Note that these histograms show the likelihood of such sys- 
tems to be zero at redshifts z ^ 10. However, this is simply because 
the Millennium simulation merger trees do not resolve progenitors 
at these redshifts. As mentioned before, we have minimized the ef- 



fect of this shortcoming on our results by requiring that the SMBHs 
always follow the M - cr relation initially. 

High mass galaxies (Mq « lO''' M©) are more likely to have 
multiple BHs in their nuclei at higher redshift. About 60% of these 
galaxies have more than 2 BHs between redshifts z~ 2 and 10. This 
fraction is less than 40% for the low mass galaxies (Mo a 10'* Mq) 
The likelihood of occurrence of more than 3 and 4 BHs is similar, 
about 30%, in the two categories of simulation. However, for the 
high mass galaxies this likelihood is spread out over a wider range 
in redshift, again due to the higher rate of major mergers. 

It is extremely rare for Milky Way-sized galaxies (halo mass 
Mo ~ 10'^ Mq) to have more than three SMBHs in their nuclei at 
any moment in their assembly history. Indeed, in our simulations 
of these galaxies, only one run shows a triple BH system. The main 
reason behind this is the smaller number of major mergers for these 
galaxies. Moreover, it is easier for SMBHs to escape the nuclei of 
predominantly small mass progenitors of these galaxies. 



6.4 EfTects on the stellar distribution 

Most bulges and early-type galaxies have a shallow cusp near their 
centre. The mass distribution in this region can be described as 
a po wer law p oc r~^. Most galaxies have slo pe 0.5 < y < 
2.0 (jFerrarese et al.ll2006l ; JMerritt & Szellll2006h . We expect the 
constituent SMBH in the bulge to affect the mass distribution 
within its radius of influence. Only two galaxies, the Milky Way 
(Genzeletal. 2003) and M32 JLauer et alj|l998[) . have been re- 
solved at these small distances. Both these galaxies have y « 1.5 in 
their innermost regions. 

It is commonly postulated that cores can form in elliptical 
galaxies and spiral bulges due to mass ejection by a hard binary 
SMBH (e.g. Milosavlievic & Merritt 2001). However, the mass 
ejected by a hard binary is of the order of the black hole mass. 
In other words, the mass deficiency M^^f, which is the difference 
between the mass of the initial and final density distribution in a re- 
gion around the centre, is roughly M\,h, the total mass of the SMBH 
binary. The possibility of enhanced mass deficit because of re- 
peate d core passages of recoiled bl ack holes ([Gualandris & MerrittI 
120081) and due to repeated mergers ( lMerritll2006l) has been consid- 
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Figure 14. Mass deficiency versus number of coalescences averaged over 
ten simulation runs. The presence of multiple SMBHs generally leads to 
larger mass deficiency compared to a single hard SMBH binary. 



ered in the literature. Our simulations allow us to understand the 
effect of both of these factors in addition to the mass deficit pro- 
duced by simultaneous presence of multiple SMBH in the galactic 
bulge. 

Figure |6^ shows the cusp evolution in two of our simulations, 
each of which has four SMBHs and three coalescences. Density 
profiles after each coalescence is shown. Strong core formation is 
clearly seen. We calculate Mjef/Mbh for ten such runs and show 
the average result in Figure. Clearly Mdcf /Mbh is much larger when 
multiple SMBHs are present. Values of Mjef /Mbh ~ 5 have been 
observed in large elliptical galaxie s tGraham.2004l ; lFerrarese et al.l 
120061 : iHopkins & HemquistlllOlCI) . Our model explains the occur- 
rence of such systems. Since the star-star relaxation time in large 
elliptical galaxies in expected to be ~ 10'" yr, we can expect them 
to carry the signature of core formation at high redshift due to mul- 
tiple SMBHs. At lower redshift our simulation are applicable to 
spiral bulges, which have much lower relaxation time scale (~ 10' 
yr). Indeed in the runs where a single black hole is left for z < 2, we 
find the formation of a Bahcall-Wolf cusp. This is consistent with 
the observed structure of the Milky Way bulge. 

The above considerations regarding cores in galaxy luminosity 
profile are also applicable to dark matter cores. The ejection of dark 
matter particles by the black holes will produce a core similar in 
size to the stellar core. 



7 OBSERVATIONAL SIGNATURES 

From the results of our simulations described above, we expect 
about 30% of the galaxies within haloes with a present-day mass 
of Mo a 10''' Mq to contain more than two SMBHs at redshifts 
2 < z < 6. For more massive haloes with Mq > 10" Mq, this frac- 
tion is almost 60%. However, since few such systems have been 
unambiguously observed so far, we consider some observational 
signatures that would indicate their existencqj. Apart from their ef- 
fect on the stellar mass distribution, multiple SMBH systems lead 

^ Some systems with triple active gala ctic nuclei (A GNs) were reported so 
far. Examples are NGC 6166 and 7720 iTonrvll984[) and SDSSJ1027+1749 



to an enhanced rate of tidal disruption of stars, modified gravita- 
tional wave signals compared to isolated BH binaries, and slingshot 
ejection of SMBHs from galaxies at high speeds^ 

From the results of scattering experiments. Ichen et alj ( 12003) 
found that the stellar tidal disruption rates due to three-body in- 
teractions between a hard, unequal-mass SMBH binary with fixed 
separation and a bound stellar cusp is higher by several orders of 
magnitude than the corresponding rates for a single SMBH. In par- 
ticular, they find that the stellar tidal disruption rate is about 1 yr"' 
for an isothermal stellar cusp with cr = 100 km s"' containing an 
SMBH binary of total mass 10^ M©. In comparison, the correspond- 
ing rate for a single 10^ Mq black hole is about 10 * yr ' . The dura- 
tion of the tidal disruption phase is about 10' yr. This enhancement 
in the tidal disruption is du e to the Kozai-Li dov effect and due to 
chaotic resonant scattering JChen et al.l20Il '). Tidal disruption of a 
star results in about half of the stellar mass being inserted in bound 
elliptical orbits. When it falls back in the black hole, this mass gives 
rise to a bright UV/X-ray emission ("tidal flare") lasting for a few 
years. One such event may have already been recently observed in 
the form of high-energy transients that can be mode l ed as sudden 
accre t ion events onto an S MBH I ILevan et alj|201 Ik iBloom et al.l 
|20I lUZauderer et al.ll201 ih . 

We expect similar enhancement in the rate of stellar tidal dis- 
ruption in systems with multiple black holes. Firstly, the presence 
of multiple SMBHs increases the combined tidal disruption cross 
section of the black holes. (Although this will only enhance the 
tidal disruption rate by a factor of a few.) Secondly, even before 
they closely interact, the presence of a third SMBH affects the tidal 
disruption event rate onto an SMBH binary by scattering stars into 
the binary's loss cone at a r ate that increases as inve rse square of its 
separation from the binary JHoffman & Loebl2007h . Thirdly, as we 
saw above, multiple SMBH systems are likely to contain recoiled 
black holes, which have been kicked either due to anisotropic grav- 
itational wave emission after coalescence, or due to the gravita- 
tional slingshot. Sudden recoil promptly fills the loss cone of these 
black holes. The resultant enhancement in the tidal d isruption event 
rate c an be substantial, increasing it up to 0.1 yr"' JStone & Loebl 
I2OIII) . Furthermore, if their recoil velocity is not too high, these 
recoiled SMBHs oscillate around the stellar core with decreasing 
amplitude due to dynamical friction. This motion results in their 
repeated passages through the stellar core, thereby increasing the 
stellar tidal disruption event rate. 

Another observational signature of systems with multiple 
SMBHs is gravitational waves (GWs). The GW emission from 
binary and triple SMBHs has been stud i ed in the literature 
('Wvithe & Loebl |2003J : ISesana et al.l |2004| : lAmaro-S eoane et al.l 
2010). Space-based detectors like the Laser Interferometer Space 
Antenna (LISA) are expected to be sensitive in the frequency 
range ~ 10"''-10"' Hz. This corresponds to the inspiral of SMBH 
systems with total mass ~ 10 '' - 10'° Mg. Puls ar timing arrays 
(PTAs) li ke the Parkes PTA JManchesteil |2008|) and the Euro- 
pean PTA jjanssen et al.l2008l) and ground-based detectors like the 
North American N anohertz Observatory for Gravitational Waves 
djenet et alj 120091) are sensitive to even lower frequencies of ~ 
IQ-s 10-6Hz. 

lYunes et alj ( 1201 ih studied modifications due to the presence 



i Liuet al.ll201ll) . The first two objects are cD galaxies at z ~ 0.03 and 
the latter is at z w 0.06. All three are kpc-scale triples. It is possible that 
NGC 6166 is simply a super position of a centr al cD galaxy and two low- 
luminosity elliptical galaxies ILauer et all 19981) . 
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of a secondary SMBH in the waveform of an extreme mass-ratio 
inspiral (EMRI) of a stellar mass objects into an SMBH. They 
find that a 10* M© SMBH will produce detectable modifications 
if it is within a few tenths of a parsec from the EMRI system, al- 
though this distance increases for higher mass SMBHs. In this pa- 
per, we have quantified the presence of such 'massive perturbers.' 
The resultant modifications to gravitational waveforms will be a 
distinct signature of multiple-SMBH systems. Futhermore, such 
systems often contain binaries that have phases of very high ec- 
centricities, created via mechanisms like the Kozai-Lidov effect 
JHoff^man & LoebluOOTn . Such binaries are expected to to emit in- 
tense b ursts of high-frequency grav itational waves at the orbital pe- 
riapsis f Amar o-Seoane et alj20 1Cr). As a result, sources that would 
normally emit outside of the frequency windows of planned grav- 
itational wave searches may be shift ed into observable range. For 
example, lAmaro-Seoane et alj ( 120101 find that a few to a hundred 
gravitational wave bursts could be produced at a detectable (1 ns) 
level within the PTA frequency range if the fraction of SMBH 
triplets is > 0. 1 . 

Presence of triple SMBHs also has important implications for 
gravitational wave searches using match ed-filtering by possibly re- 
quirin g additional waveform templates JAmaro-Seoane & Freitad 

i2onh . 

Lastly, an observable signature of these systems will 
be the presence of w andering SMBHs in the large haloes 
JHoflFman & Loebll2007h . We have shown that about 10% of the 
SMBHs are ejected at velocities > 2000 km s "' due to the sling- 
shot mechanism. This high-speed black holes will spend 1-10 Gyr 
in the outskirts of the halo. However, it is not clear whether detect- 
ing this population of wandering black holes will be possible. 



8 CONCLUSIONS 

In this work, we have addressed the formation of galactic nuclei 
with mutiple SMBHs. We performed accurate N-body simulations 
of mergers of galactic nuclei with SMBHs in a cosmological set- 
ting. Our calculation uniquely incorporated cosmological mergers 
of galaxies with an accurate treatment of dynamical interactions be- 
tween SMBHs and stars, which we achieved using the direct sum- 
mation N-body code, NBODY6. T he need for such simulations ha s 
been recognized in the literature JMerritt & Milosavlievid 120050 . 
Our main conclusions are as follows: 

• In the absence of gas, high mass galaxies (Mo > 10'* M© 
at z = 0) are generically expected to have had multiple SMBHs 
in their nuclei during their assembly history. Our simulations sug- 
gest that ~ 30% galaxies within haloes with a present-day mass of 
Mo X lO''* Mo (M-6 a 10" Mo) contain more than two SMBHs at 
redshifts 2 < z < 6. For more massive haloes, with Mo > 10'^ Mo 
(M,=6 ~ 10'^ Mo), this fraction is almost 60%. This is in contrast 
to lower-mass galaxies (Mo « 10'^ Mq; M_6 x 10'" Mq), which 
rarely host more than two SMBHs in their nuclei at any moment in 
their assembly history. 

• High mass galaxies as well as their low mass counterparts 
are rarely expected to retain more than two SMBHs in their nuclei 
at the present epoch. SMBH coalescence and ejection reduces the 
number of SMBHs on the time scale of a Gyr. Furthermore, major 
mergers are rare at lower redshift. We also find that the number of 
SMBHs in galactic nuclei is rarely reduced to zero at z = 0. Less 
than 5% of our high-mass runs resulted in such galaxies. 

• SMBH coalescence is common at high redshifts. Subsequent 
recoil due to anisotropic gravitational wave emission often results 



in escaping SMBHs. Some of these SMBHs add to the wander- 
ing population of black holes in the galactic halo. In a few cases, 
this process also results in galactic nuclei with no SMBH near their 
centres. BH-BH interaction also leads to ejected SMBHs via the 
slingshot mechanism. While most of ejected SMBHs have veloc- 
ities < 500 km s"', about 10% SMBHs are ejected at very high 
velocities exceeding 2000 km s"' . We also find binary SMBH ejec- 
tion in < 10% of the cases. 

• Multiple SMBHs have a strong effect on the stellar distribu- 
tion due to three-body interactions and core passages. The resulting 
mass deficit is usually much larger than that due to a single SMBH 
binary because of resonant BH-BH interactions and GW recoil of 
the BH remnant. We observe long-term oscillations of the BH-core 
system that could explain observations of offset AGNs. Th is has 
implications for recent observations by Icivano et al.l ( 120101) of a 
z = 0.359 system that potentially contains a recoiled BH. 

• The presence of multiple SMBHs will have important effects 
on the rate of tidal disruption of stars in galactic nuclei due to en- 
hanced tidal disruption cross section, scattering of stars by other 
BHs, prompt loss cone refilling due to GW recoil and gravitational 
slingshot. Similarly, the presence of more than two BHs in a hier- 
archical triple is expected to leave a signature in the GW emission 
from the inner binary. This signature could be observed with future 
GW observatories, such as LISA. Finally, we also expect such sys- 
tems to give rise to a distinct population of wandering SMBHs that 
could travel in large haloes over long time scales of a few Gyrs. 

The presence of gas could alter the above picture to some ex- 
tent. However, simulations of binary BHs in gaseous environment 
have not reached sufficient resolution to confirm this. Moreover, 
we expect that at high redshifts, AGN activity triggered by galaxy 
mergers could efficiently drive gas away from the shallow potential 
wells of the galaxy. Our work can also be extended by calculating 
late stages of binary SMBH evolution more consist ently. New reg - 
ularization techniques to do this are now available ( Aarseth 2003'); 
we defer their use to future work. Furthermore, multiple SMBH 
systems can als o form in additional way s, for example by fragmen- 
tation of disks JGoodman & TanI 20041). How ever, these systems 
would evolve by migration jKocsis et al.ll201ll) on a much shorter 
time scale than considered here. 
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